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Abstract In its simplest statistical-mechanical description, 
a granular fluid can be modeled as composed of smooth in- 
elastic hard spheres (with a constant coefficient of normal 
restitution a) whose velocity distribution function obeys the 
Enskog-Boltzmann equation. The basic state of a granular 
fluid is the homogeneous cooling state, characterized by a 
homogeneous, isotropic, and stationary distribution of scaled 
velocities, F(c). The behavior of F(c) in the domain of ther- 
mal velocities (c ~ 1) can be characterized by the two first 
non-trivial coefficients («2 and 03) of an expansion in So- 
nine polynomials. The main goals of this paper are to review 
some of the previous efforts made to estimate (and measure 
in computer simulations) the a-dependence of ai and 03, 
to report new computer simulations results of 02 and 03 for 
two-dimensional systems, and to investigate the possibility 
of proposing theoretical estimates of 02 an d a i with an opti- 
mal compromise between simplicity and accuracy. 

Keywords Homogeneous cooling state • Sonine coeffi- 
cients • Linear approximations 



1 Introduction 

The prototype model for a granular gas is a system of smooth 
inelastic hard spheres characterized by a coefficient of nor- 
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mal restitution < a < 1 |Q~|[2). Kinetic theory arguments 
similar to those followed in the elastic case allow one to de- 
rive the Boltzmann and Enskog equations [ 3 1 for the velocity 
distribution function f(r, y,t ). 

Perhaps the basic and simplest physical state for a gran- 
ular gas is the homogeneous cooling state (HCS), where the 
gas is isolated and has an isotropic and spatially uniform 
distribution 0. In this state, the collisional loss of energy, 
characterized by the cooling rate £, makes the mean kinetic 
energy (directly related to the so-called granular temperature 
T) decay monotonically in time following Haff's law [4|: 



T(t) = 



no) 



[i + iC(o>] 



(i) 



Therefore, the distribution function evolves in time toward 
a delta function, i.e., /(v) — * «5(v), where n is the num- 
ber density. However, the simplicity of this trivial asymp- 
totic state is deceiving since the distribution function actu- 
ally reaches an interesting scaling (or self-similar) form 



/(V) - nv Q d (t)F(c(t)), c(0 = v/v (f). 



(2) 



Here, d is the dimensionality and vo(f) is the thermal speed 
defined by 



dTV 2 /(T,0 



By definition, 

where the (scaled) velocity moments are 
(c 2 ") = / dcc 2p F(c). 



(3) 



(4) 



(5) 
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In the HCS, the Enskog-Boltzmann equation for the prob- 
ability distribution function F(c) of the reduced velocity is 

m 

ju 2 d 



, cF(c)=I[c\F,F], 
where the collision operator is 

/[ci|F,F] = Jdc 2 j dff0(ci2-ff)(ci 2 -ff) 

x[a- 2 F(c'()F(^-F( Cl )F( Cl )] 
and we have introduced the collisional moments 

H2 P = - I dcc 2 Pl[c\F,F], 



(6) 



(7) 



(8) 



so that 2^,2 /d is the dimensionless cooling rate. In Eq. (0, 
C12 = Ci — C2 is the relative velocity of the colliding par- 
ticles, a is a unit vector directed along the line of centers 
from the sphere 1 to the sphere 2, is the Heaviside step 
function, and ((;'/, c 2 ') are the precollisional or restituting ve- 
locities yielding (ci ,02) as the postcollisional ones, i.e. 



1 



•'[.2 = t\,2T^+a )(ci2 • g)g . 



(9) 



The exact solution of Eq. (|6]l is not known. Except in the 
elastic case (a — 1), the Maxwellian F(c) — n~ d l 2 e~ c ~ 
(c) is not a solution. In particular, if a < 1 , it is known that 
F(c) develops an exponential high-energy tail Q6), 



F(c). 



(10) 



A convenient way of characterizing the deviation of F(c) 
from (j) (c) in the regime of low and intermediate speeds is 
through the Sonine polynomial expansion 



F(c) = 0(c) 



1 



k=2 



V) 



(11) 



where L\ are generalized Laguerre (or Sonine) polynomi- 
als 0. The first two non-trivial coefficients are 02 an d a 3- 
They are related to the fourth and sixth velocity moments as 



(c 4 ) 



d(d + 2) 



(l+a 2 ), 



6 = d(d + 2)(d + 4) 



(1 + 3«2 — a i) 



(12) 



(13) 



The Sonine expansion (fTTT i is known to be only asymp- 
totic [8 1, so that its practical applicability is restricted to low 
and intermediate velocities (say c < 1), in which case the 
most relevant coefficients are a-i and 03. Therefore, the de- 
termination of these two coefficients is important to quantify 
the basic deviations of the HCS distribution F(c) from the 



Maxwellian 0(c), at least for c < 1. This explains the inter- 
est this problem has attracted over the years 1 2, 5, 8, 9, 10, 1 1 , 
[Hl[l3l[T4l[T5l[T6l[T7l[T8l . Apart from this formal motivation, 
the knowledge of 02(a) and 03(a), especially in the case 
of the former, is needed to evaluate the dependence of the 
transport coefficients on inelasticity ||8l [T9ll20l . 

The aim of this paper is three-fold. First, some of the ef- 
forts done in the last dozen years or so to estimate 02 an d 
03 by theoretical tools and to measure them in simulations 
are briefly reviewed in Sec. [2] Next, we explore the possi- 
bility of deriving theoretical expressions for ci2 and 03 with 
an optimal compromise between simplicity and accuracy. To 
that end, we restrict ourselves to the class of linear approx- 
imations, revisit some of the ones already proposed in the 
literature, and construct a few new ones in Sec. [3] Those ap- 
proximations are compared with new (d = 2) and recently 
published [ 17 1 (d = 3) computer simulations in Sec. [4] Sec- 
tion [5] addresses the case of granular gases heated with a 
white-noise thermostat. Finally, the conclusions are summa- 
rized in Sec.|6j while some complementary material is rele- 
gated to the Appendices. 



2 A brief review of previous results 

Taking (even) moments in both sides of Eq. (|6]l one gets the 
exact set of moment equations 



-2/n 



P>2, 



(14) 



where use has been made of Eq. (0). The condition p > 2 is 
introduced because Eq. (fl4l becomes an identity for p = 
and also for p = 1 . 

It is important to notice that the collisional moments are 
functionals of the distribution function, so that Eq. dT4b im- 
plies a coupling among all the Sonine coefficients a^ and 
there is no a priori reason to expect a chosen truncation to 
provide accurate results for a subset of coefficients. On the 
other hand, most of the routes followed to get theoretically 
based results assume some kind of truncation and/or order- 
of-magnitude simplification. More specifically, one usually 
approximates the first few collisional moments jj.2 P by in- 
serting the expansion ( fTTT i into Eqs. (0 and ((8), truncating 
the expansion after a certain order and, in some cases, ne- 
glecting nonlinear terms. The resulting set of approximate 
equations is then solved algebraically to obtain estimates for 
the desired coefficients a^. In principle, these estimates are 
uncontrolled and can be assessed only after comparison with 
computer simulations. 

The first application of this method was carried out by 
Goldshtein and Shapiro in a pioneering and extensive paper 
|9l . They derived a simple expression for «2 in the three- 
dimensional (d = 3) case by taking a linear approximation 
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(namely, neglecting a\ and with k > 3) in Eq. (TBI) with 
p — 2. Here we quote their result: 



«2 : 



16(l-a)(l-2a 2 ) 
Ao + Aia + A2a 2 (l -a) : 



(15) 



where (Ao,Ai,A 2 ) = (401,-337,190). They noted that, ac- 
cording to their estimate, the magnitude of «2 was quite 
small (|fl2| < 0.04). However, there was a small algebraic 
mistake in their derivation that was subsequently corrected 
by van Noije and Ernst [5 1, who also generalized the result to 
arbitrary d. The expression derived by van Noije and Ernst 
(vNE) maintains the form ( fT3l >. except that (Ao,Ai,A2) = 
(9 + 24c/, 8c/ — 41,30), what in the three-dimensional case 
becomes (Ao,Ai,A2) = (81,-17,30). This yields values of 
|«2 1 up to five times larger than those predicted by the (wrong) 
original expression by Goldshtein and Shapiro. Although 
published in 1998, vNE's result had been circulating ear- 
lier and so in 1996 Brey et al. iflO j confirmed its accuracy 
for d = 3 and a > 0.7 by comparison with DSMC sim- 
ulations 1271 of the Boltzmann equation. Brey et al. iflOl 
also presented simulation data for (c 6 ) but they did not ex- 
tract from them the corresponding values of 03. When this is 
done from Figs. 5 and 6 of Ref. ifTUl . one gets a$ ~ —0.005 
for 0.7 < a < 0.9 and d — 3. More recently, Ahmad and 
Puri |[T8l carried out large-scale event-driven molecular dy- 
namics simulations and measured aj_ and 03 in the HCS for 
a > 0.7 in two- and three-dimensional systems. The results 
confirmed the accuracy of vNE's expression of aj_ in that 
range of inelasticity and showed that 03 was typically four 
to five times smaller than 02. These authors also studied the 
time evolution of for k = 2-5 to monitor the transition 
from the HCS to an inhomogeneous cooling state. 

In 1999, Garzo and Dufty [11] studied the HCS for three- 
dimensional binary mixtures. Neglecting again a\ and 
with k > 3, they obtained explicit expressions for the Sonine 
coefficient 02 of both species as functions of the three coef- 
ficients of restitution and of the temperature ratio Ti/T^. To 
close the problem, it is necessary to determine T1/T2 from 
the condition of equal cooling rates for both species, yield- 
ing a highly nonlinear equation. The results showed that typ- 
ically the component made of particles with a larger mass 
has a higher temperature and a higher value of 02. The the- 
oretical predictions were later validated by DSMC simula- 
tions lfT5ll . 

The authors pointed out in 2000 lfl2ll that the linear ap- 
proximation to get d2 is not univocally defined, the result 
depending on the way that the quantities in Eq. (TBI) are ar- 
ranged. In particular, if Eq. (fT4l > for p = 2 is rewritten as 
fj,4/(c 4 ) = 2jtt2/ '(c 2 ) and then a linear approximation is ap- 
plied the result is given again by Eq. ( TT3T > but with (Ao, Ai , A2) = 
(25 +24c/, 8c/ - 57 , -2), implying (Ao , A] , A 2 ) = (97 , -33 , -2) 
for d = 3. This alternative expression is hardly distinguish- 
able from vNE's if a > 0.5 but provides up to 16% smaller 



values than the latter for higher inelasticities. We performed 
DSMC simulations of 02 for d = 3 and a > 0.2 which ex- 
hibited an excellent agreement with our alternative linear 
approximation. Moreover, DSMC data of 03 were also pre- 
sented in Ref. 03- While c? 3 ~ -0.005 for 0.6 < a < 0.9, 
a relatively rapid decay of 03 for higher inelasticities was 
observed. 

Brilliantov and Poschel |[T3l were possibly the first ones 
to depart from the linear approximation. They neglected 
with k > 3 but retained a 2 in Eq. (TT~4T > with p = 2, thus ob- 
taining a closed cubic equation for «2 ( an d d — 3). Its phys- 
ical root is very close to vNE's expression for a > 0.4 but 
becomes up to 10% larger than it for smaller values of a. 
Taking into account the simulation results presented in Ref. 
fl2l . it turns out that the physical root of the cubic equation 
deviates in the wrong direction from vNE's simpler linear 
approximation. This paradoxical outcome shows that the So- 
nine coefficients with k > 3 are not negligible for a < 0.4. 

A different truncation scheme was followed by Huth- 
mann et al. [ 14|, who assumed that — G{e k ), where £ ~ 
\a2^ 2 was treated as a small parameter. The solution to or- 
der k > 2 was obtained by taking Eq. (TPfl i for p = 2, . . . ,k 
and formally neglecting terms of order e with i > k. The 
second-order solution recovers the vNE result for 02. In the 
third-order solution one has a set of two linear equations 
for «2 an d fl 3> but the fourth-order approximation involves 
a set of three equations for a%, 03, and 04 that include a 2 , 
so that the problem becomes nonlinear for k > 4. In this 
approach, the coefficient is renormalized as the trunca- 
tion order increases. Huthmann et al. applied their scheme 
to d = 2 and observed that the values of «2 obtained from 
order e 2 to order e 6 were practically the same as long as 
a > 0.6. However, those values were dramatically sensitive 
to the truncation order for higher inelasticities, thus indicat- 
ing that the assumption = ff{e k ) fails if a < 0.6. Molec- 
ular dynamics simulations showed a good performance of 
vNE's expression for hard disks and a > 0.4. 

Coppex et al. ITBI tried an approach to estimate 02 dif- 
ferent from those based on Eq. ([Pil l. They started from Eq. 
(O in the limit c — > and then inserted the Sonine expan- 
sion (fTTT > by neglecting a 2 and with k > 3. The solution 
of the resulting linear equation for ci2 had the structure of 
a polynomial of fourth degree in a 2 divided by a polyno- 
mial of eighth degree in a (with no a 5 and a 7 terms). Al- 
though promising, this alternative method yields poor results 
for small and moderate inelasticities and only improves over 
the vNE benchmark formula if a < 0.4, as comparison with 
DSMC data for d = 2 shows QH. Coppex et al. also elab- 
orated further on the ambiguity of the linear approximation 
for «2 pointed out in Ref. Ifl2l . In particular, they showed 
that the linear approximation as applied to [14/ (c 4 ) = 2\i%j (c 2 ) 
and to \li,{c 2 ) /2/X2(c 4 ) = 1 presented a very good agreement 
with their two-dimensional simulations. 
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More recently, Brilliantov and Poschel 1171 have con- 
sidered the linear approximation of Eq. ([Pil l with p = 2 and 
p = 3 by neglecting a\, «2«3, flf, and with > 4. This 
gives a set of two linear equations for a 2 and 03 that is ac- 
tually equivalent to Huthmann et al.'s method to order e 3 . 
By comparing with their own DSMC simulations for d = 3, 
Brilliantov and Poschel showed that their expression of a%, 
while rather more complicated than vNE's, provided worse 
estimates than the latter for a < 0.6. As for their expression 
of a 3 , it was quite good for a > 0.6 and exhibited the right 
qualitative behavior for larger inelasticities. Apart from a 2 
and «3, they also measured 04, as, and a^ in the simulations, 
showing that their values were not negligible if a < 0.6. In 
fact, these authors argued that the Sonine expansion breaks 
down for large inelasticity. 

Using the asymptotic high-velocity tail ( TTOl ). Noskowicz 
et al. [8] have shown that a k oc (-A/E, 2 ) k (k + 1)! for large 
k, so that the series ( fTTT i is divergent, although it is asymp- 
totic and Borel resummable. On the other hand, the Sonine 
expansion of the modified function Gy(c) = e+'^T) 6 ' F(c) 
does converge for < 7 < 5. Truncating the Sonine ex- 
pansion of G r (c) after k = N s (with N s = 10, 20, and 40), 
Noskowicz et al. obtained numerically the associated So- 
nine coefficients a k ", k = 0, 1, . . .N s , with the help of sym- 
bolic software. Once Gy(c) is (approximately) determined in 
this way, the Sonine coefficients a^ of the true distribution 
function F(c) = e( 1- ^ 6 ' Gy(c) can be obtained by quadra- 
tures. The numerical results for a 2 , which are well fitted in 
the three-dimensional case by Eq. (Q3) with (Ao,Ai,A2) = 
( 104. 1 , —5 1 .43 , 78 .67), confirmed that vNE's expression over- 
estimates a 2 for a < 0.5, while the alternative expression 
proposed in Ref. |[T2l is rather accurate (although it slightly 
overestimates a 2 ). 



3 Theoretical estimates from linear approximations 

Our main goal now is to get estimates of the Sonine coef- 
ficients a 2 and a 3 by the application of approximations that 
neglect the coefficients a# with k > 4 as well as the nonlinear 
terms a\, a 2 a 3 , an d ^3- As we will see, this recipe is not a 
systematic method and so it does not provide a unique result. 

Given a functional X [F] of the scaled probability distri- 
bution function F(c), henceforth we will use the notation 
^a 2 ,a 3 {X} t0 denote an approximation to X[F] obtained by 
using Eq. ( fTTl i and neglecting a^ with k > 4 and nonlinear 
terms (like a\, a 2 a 3 , and af). Furthermore, if 03 is also ne- 
glected, the corresponding approximation will be denoted 
by Jz? fl2 {X}. In particular, in the case of the collisional mo- 
ments defined by Eq. ^ with p = 1,2, and 3, one gets 



Xi 2 ,a 3 {^2} = A Q +A 2 a 2 +A 3 a 3 , 



■^a 2 ,a 3 {^4} =-60 + ^202+^303, 



(16) 



(17) 



^a 2 ,a 3 {/J-6} — Co + C 2 a 2 + C 3 a 3 ■ 



(18) 



The expressions for the coefficients A,, Bj, and C, as func- 
tions of a and d were derived by van Noije and Ernst 
and by Brilliantov and Poschel lT2T fT7l . They are given in Ap- 
pendix [A] Obviously, «5f a , {fl 2 } and 5£ Ul {/X4} are obtained 
by formally setting A3 — > and B 3 — > on the right-hand 
sides of Eqs. ( [TBI ) and ( fT71 ). respectively. 

The exact equation ( fT4b becomes an approximation when 
it is linearized with respect to a 2 and 03. For p = 2 and p = 3 
we get 



= X n ,a 3 <j M4 -2^2 -73T 



B - (</ + 2)A + [B 2 - (a! + 2) (A + A 2 )] a 2 
+ [B 3 -(d + 2)A 3 ]a 3 , 



(19) 



= « 3 |/^6 — 3jU2-^3y 

= C -^(d + 2)(d + 4)A +[C 2 -~(d + 2)(d + 4) 

x (3A + A 2 )]a 2 + [C 3 -^{d + 2) (d + 4) (A 3 - A )]a 3 . 

(20) 

The non-systematic character of the linearization method 
is made evident if one proceeds in the same way, except that 
Eq. (O is rewritten in the equivalent form 



M2 . 



(21) 



This equation shows that the rescaled collisional moment 
ji 2p / (c 2p ) is just proportional to p. Now, instead of Eqs. JT91 
and ( f20b we have 

U - 1 ( C 4) / (c 2 ) 

= B -(d + 2)A + [B 2 -B Q -(d + 2)A 2 ]a 2 

+ [B 3 -(d + 2)A 3 ]a 3 , (22) 



= JS?* 



M6_ 



M2 



= C --(</ + 2) (</ + 4)A + [C 2 -3C - j(d + 2) 



x(d + 4)A 2 



a 2 - 



C 3 +C Q ~-{d + 2)(d + 4)A 3 



a 3 . 
(23) 

Note that Eq. (1221 is obtained from Eq. $1% if one formally 
replaces (d + 2)Aq by B$ in the coefficient of a 2 . Likewise, 
Eq. ( 1231 is obtained from Eq. ( f20T > by formally replacing 
j{d + 2){d + 4)Ao by Co in the coefficients of a 2 and 03. 
Nevertheless, the approximations (Tf9] l and ( f20b are differ- 
ent from the approximations ( 122b and ( 1231 . respectively, so 
they provide different estimates of the coefficients a 2 and 
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a?,. Henceforth we will refer to the approximations ( fT9b and 
d20l i. which are based on Eq. (fT4l i. with the label "a". Analo- 
gously, the approximations (f22b and d23l l. which are based 
on Eq. d2Tl i. will be labeled by "b". Of course, other al- 
ternative ways of rewriting Eq. ([Pil l are possible I112II16L 
With independence of the adopted approach (say "a" or "b"), 
there are two basic classes of approximations: either a 3 is 
neglected versus a 2 in the equation for /I4 (but not in the 
equation for or both Sonine coefficients are treated on 
the same footing. 

3.1 Class-I approximations: 03 -C a 2 

Let us first assume that 03 can be neglected versus 02 in 
either Eq. $1% (approach "a") or Eq. d22b (approach "b"): 





B -(d + 2)A + [B 2 -(d + 2) (A +A 2 )] a 2 , (24) 



- ^"2 1 / 4\ L (c 2\ 



= B -(d + 2)A Q + [B 2 -B -(d + 2)A 2 ]a 2 . (25) 
These are linear equations for a 2 whose respective solutions 



are 



4» = 



4 b («) 



B Q -(d + 2)A 
{d + 2){A 2 +A Q )-B 2 

16(1 -a)(l -2a 2 ) 
9 + 24c/- (41 - 8c/)a + 30(1 - a)a 2 ' 

B Q -{d + 2)A 
(d + 2)A 2 -{B 2 -B Q ) 

16(1 -a)(l -2a 2 ) 



(26) 



(27) 



25 + 24d- (57 -8c/)a- 2(1 -a)a 2 ' 
where in the last steps use has been made of the explicit 
expressions of Ao, A 2 , Bo, and B 2 . As recalled in Sec. [2] 
the method labeled here as "la" was the one first followed 
by Goldshtein and Shapiro [9 |, the corresponding estimate, 
Eq. (1261 . being first obtained by van Noije and Ernst [5 1. The 
alternative method "lb", Eq. d27i>. was proposed in Ref. [12|. 
It is interesting to note that 



4 b = 



(28) 



Comparison with DSMC data shows that the estimate a 2 b is 
superior to a 2 a f° r a ^ 0-6 lfl2ll . Other class-I approxima- 
tions for a 2 were considered by Coppex et al. ifTBI and are 
more generally described in AppendixlBl 

Once 03 has been neglected in Eqs. ( fT9l and (1221 1. we 
can use Eq. (f20b (approach "a") or Eq. ( l23l (approach "b") 
to express 03 in terms of a 2 . The respective results are 



a 3 a (a) = G„(a,a 2 a (a)), 



(29) 



a 3 (a) 
where 

G„(a,a 2 ) 



G,(a,fl2 b («)), 



(30) 



{Co - \ (d + 2) (d + 4)A + [C 2 - 3 - (d + 2) 
x (d + 4) (3A + A 2 )j a 2 } / (d + 2) (d + 4) 
x(A 3 -A )-C 3 



(31) 



G b (a,a 2 ) = {Co-^{d + 2)(d + 4)A c 



C 2 — 3Co 



(d + 2)(d + 4-)A 2 y 2 }/\ i -(d + 2) 



x(d + 4)A 3 -C 3 -C 



(32) 



It is also possible to construct a hybrid approximation 
"Ih" in which a 2 is obtained from Eq. dZSb and 03 is subse- 
quently obtained from Eq. ( f2Qb . In that case, a 2 h = a l 2 while 



fl 3 h (a) = G a (a,4 b (a)). 



(33) 



The other hybrid possibility 03 = G/,(a,a 2 a (a)) turns out to 
be rather poor and will not be further considered here. 



3.2 Class-II approximations: 03 ~ a 2 

If «3 is formally treated as being of the same order as a 2 , 
Eqs. ( fT9l and ( f20b become a linear set of two coupled equa- 
tions for a 2 and 03 (approach "a"). This was the method re- 
cently considered by Brilliantov and Poschel fT71 . Now the 
problem is algebraically more involved than in the class-I 
approximation. The solution for a 2 is 



a 2 Ia (a) = 



where 



Ng{a) 
Da(a) ' 



(34) 



N a {a) = C 3 Bo-CoB 3 + {d + 2){A 3 Co-AoC 3 ) + -{d + 2) 
x(d + 4) [A B 3 - (A3 - A )B - (d + 2)A 2 ] , (35) 



Dfl(a) = C2B3 - C3B2 + {d+2) [(A 2 + A )C 3 - A3C2] 

+\ (d + 2) {d + 4) [(A 3 - A )B 2 - (A 2 + 3A )B 3 



+ (c/ + 2)(Ao+A 2 + 2A 3 )Ao] 
The corresponding result for 03 is 

a 3 Ia (a) = G a (a,4 Ia (a)). 



(36) 



(37) 



Note that, although the same functional form G a (cc,a 2 ) ap- 
pears in Eqs. and ([37]), obviously a^ Ia (a) ^ a!f(a). 
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The same class-II approximation can be applied to Eqs. 
(1221 and (l23T > (approach "b"). The solution is now 



af{a) 



N h (a) 
D h {a)' 



af{a) = G b {a,af{a)), 
where 



(38) 



(39) 



N b (a) = {Co + C 3 )Bo-C B 3 + {d + 2)[A 3 Co-(Co + C 3 ) 
x A ] + ^-{d + 2) (d + 4) (A B 3 - A3B0) , (40) 

D b (a) = (C2-3C )B3-(Co + C 3 )(B2-Bo) + (rf + 2) 
x[A 2 (C + C 3 )-A 3 (C 2 -3Co)] 

+l(d + 2)(d + 4)\A 3 (B 2 -Bo)-A 2 B 3 ]. (41) 

The three class-I and two class-II approximations de- 
scribed in this Section are summarized in Table Q] As said 
before, a 2 a and a 2 b = were already proposed in Refs. Q 
and Ifl2l . respectively, while a 2 a ar, d a? 3 were derived in 
Ref . IfTTll . All the other possibilities, to the best of our knowl- 
edge, have not been considered before. 



4 Comparison with computer simulations 

In order to assess the reliability of the linear estimates for 
the Sonine coefficients a 2 and 03 introduced in Sec. [3] it is 
necessary to compare them against computer simulations. 
To that end, we have performed new DSMC simulations for 
inelastic hard disks (d = 2). In the case of inelastic hard 
spheres (d = 3) we have used the extensive DSMC simu- 
lations presented in Ref. [17]. Figure [T] compares the sim- 
ulation data of a 2 with the theoretical estimates a 2 a , a l 2 h = 
a l 2 , a\ la , anc ^ a 2 b - I*- ' s observed that the best global agree- 
ment with computer simulations is provided by the two ap- 
proaches "b", i.e., the ones based on linearization of Eq. 
(|2TI) . in contrast to the two approaches "a", which are based 
on linearization of Eq. (TPfl) . Given that the class-I estimate 

a l 2 = (where 03 is neglected) is much simpler than the 8jJ.4 = /I4 — Jff a2 {^4} 1 
class-II estimate a 2 b (where 03 is retained), the former is 
preferable to the latter. In the region 0.6 < a < 1 the four 
approximations practically coincide among themselves and 
with the simulation results, thus showing that a\ and with 
k > 3 are indeed negligible in that region. On the other hand, 
our simulation data for hard disks (d = 2) show a slight im- 
provement of the two class-II approximations with respect 
to the two class-I approximations, what indicates that the in- 
fluence of a 2 is even smaller than that of 03 for 0.6 < a < 1 . 
The opposite behavior appears in the case of hard spheres 
(d = 3), although a certain scatter of the data in this case pre- 
vents us from confirming or refuting the behavior observed 
in the two-dimensional case. 



Next, we consider the third Sonine coefficient 03. The 
simulation data are compared with a l f, a lb , a lb , a 1 ^, and a^ b 
in Fig. 12 It is apparent that both approaches a lb and a 1 ^ have 
a very poor global performance, failing to account for the 
rapid increase of the magnitude of 03 when a < 0.6. How- 
ever, a good semi-quantitative agreement with computer sim- 
ulations is found for a 1 *, a^, and a" 3 . All of this implies that 
the linearization of Eq. (f2TT > with p = 3 is much less accu- 
rate than the linearization of Eq. ( TT4T > with p = 3, in contrast 
to the situation with p = 2. Interestingly enough, among the 
three estimates of 03 based on the linearization of Eq. ( TBI 
with p = 3, the best behavior is presented by the two class-I 
approximations, namely a l f for d = 2 and a lb for d = 3. As 
for the region 0.6 < 05 < 1, the two class-II approximations 
are the most accurate ones. This can be understood as a con- 
sequence of the better behavior of a 2 a and a 2 b over a l f and 
a 2 b in that region, as discussed above in connection with Fig. 
ffl 

The performance of the five linear approximations is suc- 
cintly summarized in Table Q] The best global agreement 
with simulations is achieved by the approximation "Ih", i.e., 
a 2 is autonomously obtained by linearizing Eq. (I2TI 1 with 
p = 2 and neglecting 03. Next, 03 is obtained in terms of 
a 2 and a from the linearization of Eq. ( fT4b with p = 3. The 
second best combination is "Ha", where 02 and 123 are simul- 
taneously derived from linearization of Eq. ( TT4T > with p = 2 
and p = 3. While "Ih" is simpler than "Ha", it provides a 
better estimate of 02 and 03 for high inelasticity (a < 0.6). 
This comes from a fortunate cancellation of errors and is yet 
another indication on the non-negligible character of non- 
linear terms and higher-order Sonine coefficients in that re- 
gion [17]. On the other hand, the best estimate of aj for 
0.6 < a < 1 is provided by a" 2 . 
Let us now define the deviations 



5^2 = ^2(1+02)-^ Wl +CI2)}, 



i"4 



8fi 4 = 

l+a 2 
Consequently, 



J% 2 <^ M4-2jU2- 



M4 

1 +a 2 



{d + 2)8pl 2 ~8pl4, 



(42) 



(43) 



(44) 



(45) 



(46) 



X,A^-2^ 



-2 1 , 4) 



d(d + 2) 



[{d + 2)5^2- Sju 4 ]. (47) 



The fact that a\ h = « 2 h exhibits a better agreement with simu- 
lations than ai? in the high-inelasticity region obviously im- 
plies that (d + 2)8ji 2 — 8jl4 ~ is a better approximation 
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Table 1 Summary of the linear approximations considered in this paper 



Label 



References 



Behavior of a 2 



Behavior of «3 



Equations 



"2 



"3 



0<a<0.6 0.6<a<l 



0<a<0.6 0.6<a<l 



la 


if fl2 {^ 4 -2/i2<c 4 )/<c 2 )}=0 

J? a2 ,a 3 {Ll6-3fl2(c 6 }/(c 2 }} =0 


|5| 


New 




Fair 


Good 


Good 


Fair 


Ha 


SC aim {ii4-2ti2{c A )/(c 2 )} =0 
^ 2 , a3 {Ai6-3/i 2 <t' 6 )/(c 2 )} =0 


E) 


El 




Fair 


Good 


Fair 


Good 


lb 


^{M4/<c 4 )-2/i 2 /(c 2 )}=0 
W(c 6 )-3/W<c 2 )}=0 


El 


New 




Good 


Good 


Poor 


Poor 


lib 


^,,.« 3 W<C 4 }-2^ 2 /( C 2 )}=0 
^, 2 ,a 3 {M6/<C 6 )-3M2/<C 2 )}=0 


New 


New 




Good 


Good 


Poor 


Good 


lh 


if a2 {Ai4/( C 4 )-2Ai2/<c 2 )}=0 
if fl2 , fl3 {M6-3/i 2 (c 6 >/(c 2 }}=0 


El 


New 




Good 


Good 


Good 


Fair 


Table 2 DSMC values El of a 2 , fi2, IM, 8^12, 8fi2, 8^4, and 8)14 for d 


= 3 








a 


a 2 \l2 ^4 fyi 2 




8fa 




SVa 








0.8 


-0.0141 0.8950 4.414 -0.005 




0.005 


-0.01 


-0.01 








0.6 


0.0207 1.6101 8.213 0.000 


0.000 


0.03 


0.02 








0.4 


0.0760 2.1354 11.494 0.000 


0.002 


0.09 


0.02 








0.2 


0.1274 2.4625 13.881 -0.001 


0.006 


0.20 


0.02 









than (d + 2)8Jl2 — 8/J.4 « in that region. In principle, this 
does not necessarily mean that 8jX2 ~ and 8J14 » are 
better approximations than 8J12 ~ and 8/J.4 « 0, respec- 
tively, since a certain cancellation of terms might occur in 
the difference (d + 2)8jJ.2 — 8J14. To clarify this point, Ta- 
ble|2]shows a.2, jJ-2, V-A, 8^2, 8J12, 8jX4, and 8J14 as obtained 
from DSMC simulations for inelastic hard spheres (d = 3) 
fl2l . We can observe that one typically has |5jtt2| < \8fl2\ 
and |5£Lt| < \8m\. Moreover, |5jtt2| and |5jU2| are consider- 
ably smaller than \8jM\ and | 8]m\ , thus implying that the er- 
rors made when linearizing ji4 or 1X4/ { \ +02) are generally 
larger than those made when linearizing /I2 or ^(1 +02)- 
Therefore, the property |5/I}| < \8p^\ is sufficient to jus- 
tify that the estimate of ci2 obtained by setting 8jJ>2 — * and 
8JI4 — > in Eq. (|47| | is more accurate than the one obtained 
by setting 8JI2 — ► and 8^4 — ► in Eq. d46b . 



5 White-noise thermostat 

Thus far we have assumed granular gases in the freely cool- 
ing state. The scaled quantities in this state are fully equiv- 
alent to those of granular gases kept in a steady state by a 
Gaussian thermostat [ 12 1, i.e., by the action of a determinis- 
tic non-conservative force proportional to the particle veloc- 
ity. On the other hand, a popular way of mimicking agitated 
granular gases is by means of a stochastic force assumed to 
have the form of a Gaussian white noise II22II231 . 



A simple estimate of 02 m the case of the white-noise 
thermostat was derived by van Noije and Ernst [5 1 and shown 
to be in excellent agreement with computer simulations lfl2ll . 
However, to the best of our knowledge, an analytical expres- 
sion for 03 has not been proposed so far. The aim of this sec- 
tion is to fill this gap by applying linear approximations and 
exploiting the knowledge of the coefficients A,, B,-, and C, 
in Eqs. (Q]3-([T8]) l|2ll5]fT71l . The starting point is the moment 
hierarchy, which now reads [5, 12 1 



V-2p = P 



d + 2p-2 



2p-2\ 



P>2, 



or, equivalently. 

M2p 



(c 2 p- 2 ) 



d + 2p-2 

P j M2, 



p>2. 



(48) 



(49) 



In the class-I approximation, one takes p = 2 and linearizes 
with respect to <Z2, i.e., 



3?a 2 {jl4 



-2)M 2 } = 0. 



(50) 



This condition is independent of whether we linearize Eq. 
( |48T > or Eq. d49l ), in contrast to the free cooling case. The 
solution of Eq. ( T50b is 



la lb 



Bo 



< + 2)A 



(d + 2)A 2 -B 2 

16(l-a)(l-2a 2 ) 
73 + 56t/-3(35 + 8<f)a + 30(l -a)a 2 



(51) 
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a 

Fig. 1 (Color online) Plot of the second Sonine coefficient 02 as a func- 
tion of the coefficient of restitution a for d = 2 (top panel) and d = 3 
(bottom panel). The circles represent DSMC results (d = 2: this work; 
d = 3: Ref. 1 17 1), while the lines correspond to ai? (- -•- -), o' b = 

( ), ai, Ia ( — • — )> ar >d a 2 b ( )■ ^ ne msets magnify the region 

0.6 < a < 1 



-0.12 




0.0 0.2 0.4 0.6 0.8 1.0 



a 

Fig. 2 (Color online) Plot of the third Sonine coefficient as as a func- 
tion of the coefficient of restitution a for d = 2 (top panel) and J = 3 
(bottom panel). The circles represent DSMC results (d = 2: this work; 
d = 3: Ref. 1 17 1), while the lines correspond to a l f (- -•- -), a* (- - 

-), aj 1 (- -), flj a ( — • — ), and «3 Ib ( ). The insets magnify the 

region 0.6 < a < 1 



This is the result obtained in Ref. [5|. Once a 2 is known, 03 
is determined from 







in the approximation "la" or from 



'«2*J 1 ( c 4\ 



in the approximation "lb". The results are 

a 1 f(a) = G a (a,a 1 2 a (a)), 



(52) 



(54) 



af(a) = G b (a,a 2 a (a)), 



where 



(55) 



(53) G a (a,a 2 ) = {c -hd + 2){d + 4)A + [c 2 -^(d + 2) 
x( t / + 4)(Ao+A 2 )]fl 2 }/[^(^ + 2)(rf + 4) 



xA 3 - C 3 



(56) 
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0.08 
0.06 
0.04 
53 0.02 
0.00 
-0.02 

0.06 
0.04 
0.02 
0.00 



-0.02 



n — 1 — i — ■ — i — ■ — r 




0.0 0.2 0.4 0.6 



0.8 



1.0 



a 



Fig. 3 (Color online) Plot of the second Sonine coefficient «2 as a func- 
tion of the coefficient of restitution a in the case of the white-noise 
thermostat for d = 2 (top panel) and d = 3 (bottom panel). The circles 
in the bottom panel represent DSMC results from Ref. 1121 . while the 
lines correspond to a!," = (" "*~ a 2 3 ( — * — )> ant ^ a '' b ( )• 



G b (a,a 2 ) = {c -^(d + 2)(d + 4)A + 



C 2 -Cq 



'-{d + 2)(d + 4)A 2 y 2 }/[^(d + 2) 



0.000 



-0.004 - 



-0.008 



-0.012 



-0.002 



-0.004 




-0.006 - 



-0.008 



0.0 0.2 0.4 0.6 0.8 1.0 

a 

Fig. 4 (Color online) Plot of the third Sonine coefficient aj, as a func- 
tion of the coefficient of restitution a in the case of the white-noise 
thermostat for d = 2 (top panel) and d = 3 (bottom panel). The lines 
correspond to a l f (- -•- -), af ( ), a" 3 ( — • — ), and af> ( ). 



JBor„\ _ N b( a ) JRy( n \ —r In J&tnW 

-, a 3 (a) — Ob(a,a 2 [a)), 



D b {a)- 



(60) 



where 



x(d + 4)A 3 -C 3 



(57) 



N a (a) = C 3 B - C B 3 + (d + 2)(A 3 C - A C 3 ) + - (d + 2) 



In the class-II approximations both a 2 and a 3 are simul- 
taneously obtained from 

X, 2 ,a 3 {^-{d + 2)i^ 2 }=0 (58) 

and either Eq. d52| > (approximation "Ha") or Eq. ( T53b (ap- 
proximation "lib"). The solutions are 



4>) = m> 4» = G fl (a !fl » a (a)), 



(59) 



x(d + 4)(A B 3 -A 3 B ), 

D a {a) = C 2 B 3 - C 3 B 2 + {d + 2)(A 2 C 3 -A 3 C 2 ) 
+^ (d + 2) (d + 4) [A 3 B 2 - A 2 B 3 
+{d + 2)AoA 3 -AoB 3 ], 



N b {a)=N a (a), 



(61) 



(62) 



(63) 
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D b {a) = (C 2 -C )B 3 -C 3 B 2 + (d + 2)\A 2 C 3 -A 3 (C 2 

-C Q )] + l(d + 2)(d + 4) (A 3 B 2 - A 2 B 3 ) . (64) 



Figure |3] shows the a-dependence of a!, a 



\ a" 2 , and 

ai, 10 . The three curves are very close each other, which in- 
dicates that a\ and a* with k > 3 are indeed small enough 
to be neglected in \i\ = (d + 2)ji 2 . It is interesting to note 
that a" 2 an d a 2 h we practically identical in the region 0.6 < 
a < 1, where they are slightly more accurate (at least in the 
three-dimensional case) than a\ a = a\ b - On the other hand, 
a l 2 = a l 2 and are practically indistinguishable in the re- 
gion of small a. The theoretical predictions a\*, o 3 b , a' Ia , 
and a l j h are displayed in Fig. [4] Although there are no sim- 
ulation data to compare with, it seems plausible to conjec- 
ture that the trends observed in Fig.|2]are repeated now: the 
two class-II approximations are more accurate than the two 
class-I approximations for small dissipation, a' 1 a being pos- 
sibly better than a l ^°, while the two "b" approximations are 
rather poor at high dissipation. Note that, since a!, a = a 2 h in 
the case of the white-noise thermostat, the hybrid approx- 

= a 1 ?. A feature 



af 



imation "Ih" coincides with "la", i.e. 
that becomes apparent when comparing Figs. [TJ and |2] with 
Figs.[5]and|4] respectively, is that the magnitudes of a 2 and 
03 in the white-noise case are about twice and ten times, 
respectively, smaller than those in the freely cooling state. 
This is closely related to the fact that the overpopulation of 
the high-energy tail is much smaller in the former case than 
in the latter. More specifically, instead of Eq. ( fTOb . now we 
have 



(65) 



where | is the same quantity as in Eq. (TTOb . 
6 Conclusions 



The second and third coefficients in the Sonine polynomial 
expansion of the (scaled) velocity distribution function F(c) 
of a granular gas characterize the deviation of F(c) from 
the Maxwellian and are important, for instance, in the pre- 
cise determination of transport coefficients. While for prac- 
tical purposes an interval 0.8 < a < 1 for the coefficient of 
normal restitution is sufficient, it becomes necessary from 
a more fundamental viewpoint to consider the whole range 
0< a < 1. 

The Sonine coefficients a 2 (a) and 03(a) can be mea- 
sured in computer simulations (e.g., DSMC), so one could in 
principle make a least-square fit to certain functional forms. 
However, this would not be satisfactory from a fundamental 
point of view and would not provide any insight into the in- 
tricacies of F(c) and its Sonine expansion. It is then much 
more challenging to devise theoretical approximations that 



s 



0.2 



0.1 



0.0 



-0.1 



0.3 



0.2 



0.1 



0.0 



-0.1 




V p o 




Po 



-0.2 1 ■ 1 ■ 1 ■ >— 

0.0 0.2 0.4 0.6 



0.8 



1.0 



a 



Fig. 5 (Color online) Plot of _Sf a , {^4 - (d+2)/i 2 (l +a 2 )} (o), 
^'a 2 ,aA^-(d + 2)^ 2 (l+a 2 )} (.), J? a2 {[l 4 /(l+a 2 )-{d + 2)n 2 } 
(V), and if fl2 , fl3 +a 2 ) - {d + 2)^i 2 } (T) for d = 2 (top panel) 

and d = 3 (bottom panel). The symbols are obtained from DSMC re- 
sults of a 2 and aj (d = 2: this work; d = 3: Ref. 1 17 1). The lines are 
guides to the eye. 



can be subsequently assessed by comparison with simula- 
tion results. 

In this paper we have been mainly concerned with a fam- 
ily of linear approximations to estimate 02(a) and 03 (a) in 
the freely cooling state. We have found that a good com- 
promise between accuracy and simplicity is represented by 
the hybrid approximation here denoted as "Ih", in which the 
second and third Sonine coefficients are given by Eqs. (f2Tb 
and d33b . respectively. For the benefit of the reader here we 
quote the final and complete expressions: 



02(a) 



16(l-a)(l-2a 2 



25 + 24d - (57 - 8d) a - 2( 1 - a ) a 2 



(66) 
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03(a) • 
where 



1602(a) fHCs(a) 
l-2a 2 e H cs(a) : 



(67) 



Pncs(a) = 167 + 50d- (191 + 26c/)a- 2(307 + 100c/) a 2 
+2(339 + 68c/)a 3 + 32(16 + 7d)a 4 



-32(18 + 5c/)a 5 + 144(1 - a)a 6 



GHCs(a) = 521 + 1 396c/ + 368c/ 2 - (1481 + 820cf 



-16c/ 2 )a + 4(583 + 262c/)a 2 
-20(155 + 14c/)a 3 + 280(l - a)a 4 



(68) 



(69) 



However, if more precise values in the domain 0.6 < a < 
1 are really needed, it might be preferable to consider the 
more complicated approximation "Ha" given by Eqs. d34i > 
and EJ E). 

On the other hand, It is known that in the high-inelasticity 
region a < 0.6 (i.e., once 02 becomes positive) the higher- 
order Sonine coefficients are no longer negligible ll8l fT4l[T7l . 
so that the linear approximations based on the neglect of 
nonlinear terms and of with k > 3 or k > 4 are not a 
priori reliable. This is made evident by the lack of self- 
consistency of different linear approximations used to es- 
timate C12 from the first non-trivial equation of the moment 
hierarchy, as shown in Fig. Q] What is indeed surprising is 
that the simple linear approximation d66l l provides such an 
excellent estimate both for d = 2 and d = 3. This means 
that, even though a\, 03, 04, ... are not negligible at all if 
a < 0.6, somehow they practically cancel out in the combi- 
nation JU4/ (c 4 ) — 2^2/ (c 2 ), while still playing a significant 
role in the combination m — 2ji2(c 4 ) / (c 2 ). This is clearly 
seen in Fig. [5] where we plot Jz? U2 {^4 — (d + 2)jj.2(l +02)}, 
%a 2 . n {^4 - (</ + 2)jU 2 (l +a 2 )},5?a 2 W(l +a 2 ) - (d + 2)fi 2 
and -S?a 2 ,a 3 {m/{i +^2) — (^ + 2)^2} by using the simula- 
tion data of «2 and 03. While the magnitude of the other 
three quantities rapidly increases with increasing inelasticity 
if a < 0.6, that of JSf fl2 {jU4/(l +^2) - (^ + 2)^2} remains 
small and practically constant. This property might be use- 
ful to contribute to a better understanding of the full velocity 
distribution function F(c). 

Although this paper has focused on the homogeneous 
cooling state, the analysis has been straightforwardly ex- 
tended in Sec. to a granular gas heated by a white-noise 
thermostat. In that case, the optimal combination of esti- 
mates is provided by Eqs. (fSTjl and d54l ). More specifically, 



where 

AvN(a) = 67 + 10c/-7(13-2c/)a-2(119 + 20c/)a 2 
+2(151 - 12c/)a 3 + 32(8 + 3c/)a 4 
- 32(10 + c/)a 5 + 80(l -a)a 6 , (72) 

2wn (a) = 2569 + 2932c/+ 624c/ 2 - (3529 + 2356c/ 
+240c/ 2 )a + 4(583 + 262c/)a 2 
-20(155 + 14c/)a 3 + 280(1 -a)a 4 . (73) 

It would be interesting to test the accuracy of Eq. ( T7TT > against 
computer simulations. 

A Expressions for A,, B„ and C, 

The explicit expressions of the coefficients A;, B,, and C, as functions 
of d and a are I2I5I17I 

A =^(l-a 2 ), A 2 = ^(l-a 2 ), A 3 = |-(l-a 2 ), (74) 



B =K(l-a z )[d+- + a z ), 



B 2 =K(l + a) 
K 



32 



(l-a)(l(W + 39+10cr) + (c/- 1) 



(75) 



(76) 



128 



(1 + a) [(l-a)(97 + 10a 2 ) + 2(c/-l)(21-5a)] , (77) 



Co = f (l-« 2 ) 



19 

(d + a 2 )(5 + 2a 2 ) + d 2 + — 



(78) 



3K 3 
C 2 = — (1-a 2 ) [1289 + 4(c/ + a 2 )(311 + 70a 2 ) + 172c/ 2 ] + jp, 



(79) 



}. 

c 3 = " -TT^I ( 1 " « 2 ) [ 2537 + 4 (d + a 2 ) (583 + 70a 2 ) + 236c/ 2 ] - --jS , 
1U24 Id 



3K 



_9_ 

16 
(80) 



where 



K = t= . . . , p=K(l + a)[(d- a)(3 + 4a 2 ) + 2(d 2 - a) 



y/2r(d/2) ' 



(81) 



B Other class-I linear approximations for 02 

As a generalization of Eqs. i24\ and j25) , let us consider the family of 
approximations 



CI2 



16(l-a)(l-2a 2 



73 + 56c/ - 3(35 + 8c/)a + 30(1 - a)a 2 



(70) 



(82) 



03(a) 



1602(a) AvN(a) 
l-2a 2 gwN(a)' 



(71) 



and let us denote by a* the associated solution. In particular, a 2 0,0 ' ' : 



and a. 



(1,0,0) 



: a 2 . The eight possibilities considered by Coppex et 
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al. 1161 correspond tox = 0, 1, y = 0, 1, and z = 0, 1. It is easy to check 
that the solution to Eq. 482l > is 

= (83) 

2 i+v' 

where 

^^x + y^+zf 2 -. (84) 

Equation j83t is a generalization of Eq. (28). 

Of course, other alternative possibilities exist. For instance, one 
can generalize Eq. )82t to 

where <P(X) is an arbitrary function. The corresponding approximation 
for 02 will depend on the choice of 3>(X), apart from the values of 
x,y,z. 
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